close all
clear all

system('./teta')
v=importdata('err_m.txt','\t');
b_err = v(:,1);
b_m = v(:,2);
cut_here = input('number of cols: ');
be = buffer(b_err,cut_here);
bm = buffer(b_m,cut_here);
figure
subplot 121
image(be')
title('Error')
subplot 122
image(bm')
title('Valor de m=2^k')
max(max(be))
max(max(bm))

n=size(bm);
e=reshape(be,[n(1)*n(2),1]);
m=reshape(bm,[n(1)*n(2),1]);

He=hist(e,[0:511]);
% tot1=sum(He+10^-10);
% pi=(He+10^-10)/tot1;
%     %q=beta(k)/(1+beta(k));
%     %qi=q.^[0:511]*(1-q); %Distribucion geometrica de maxima verosimilitud
%     %ti=; %Distribucion asociada al Golomb utilizado
%     %DKL(k)=sum(pi.*log(pi./ti)/log(2));
%     Hp1=sum(-pi.*log(pi)/log(2));
%     li=floor([0:511]/2^k)+k+1;
%     Lmedio(k)=sum(pi.*li);
%     MaxGain(k)=tot(k)*(Hp(k)-Lmedio(k))/(n(1)*n(2));
% figure
%     hold on
%     titl('Codificación vs. distribución empírica total')
%     plot([0:511],pi)
%     plot([0:511],2.^-li,'r')
%     legend('Distribución empírica de probabilidades','Distribución óptima asociada a la codificación utilizada')
%     hold off
for k=1:6
    eb(k,:)=e - (e+1).*(m~=2^k); %Elementos que no corresponden a este valor de k se van a -1
    Hb(k,:)=hist(eb(k,:),[-1:511]);
    beta(k)=mean(eb(k,:))^-1; %MLE de beta! 
    %Aca habria que calcular la distancia KL entre la distribucion de error (empirica) y la
    %codificacion de Golomb asociada a cada k. Podria verse si con beta se
    %podria mejorar, pero no es lo mas relevante
    tot(k)=sum(Hb(k,2:end)+10^-10);
    pi=(Hb(k,2:end)+10^-10)/tot(k);
    %q=beta(k)/(1+beta(k));
    %qi=q.^[0:511]*(1-q); %Distribucion geometrica de maxima verosimilitud
    %ti=; %Distribucion asociada al Golomb utilizado
    %DKL(k)=sum(pi.*log(pi./ti)/log(2));
    Hp(k)=sum(-pi.*log(pi)/log(2));
    li=floor([0:511]/2^k)+k+1;
    Ltest(k)=sum(2.^-li);
    Lmedio(k)=sum(pi.*li);
    MaxGain(k)=tot(k)*(Hp(k)-Lmedio(k))/(n(1)*n(2));
    figure
    hold on
    str=sprintf('Codificación vs. distribución empírica para m= %d',2^k);
    title(str)
    plot([0:511],pi)
    plot([0:511],2.^-li,'r')
    legend('Distribución empírica de probabilidades','Distribución óptima asociada a la codificación utilizada')
    hold off
    %dict(k) = huffmandict([0:511],pi); %Genero codificacion de Huffman para comparar
end

disp('Ganancia maxima (bits/pixel)')
sum(MaxGain)



figure
subplot 211
hold on
hist(e,[0:1:512])
lambda=mean(e)^-1;
aux=[0:1:512];
plot(aux,lambda*n(1)*n(2)*exp(-lambda*aux))
title('Histograma de los errores de prediccion')
hold off
subplot 212
hist(log(m)/log(2),8)
title('Histograma del parametro k')

figure
subplot 211
hold on
e2=e.*(m<128);
hist(e2,[0:1:512])
lambda2=mean(e2)^-1;
aux=[0:1:128];
plot(aux,lambda2*sum(sum(m<128))*exp(-lambda2*aux))
title('Histograma de los errores de prediccion')
hold off
subplot 212
hold on
e3=e.*(m==128);
hist(e3,[0:1:512])
hold off



